!--------1---------2---------3---------4---------5---------6---------7---------8
module material_parameters 
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind
  use global_parameters, only : irmax
  implicit none
  private
  public :: material_init, material_check, mn, ir, iw, ice_water_rock, radius_rock &
    & , iw_old, ir_old, melt, rock_water, dr_rock &
    & , radius_water, ii, radius_ice, drock, ii_old, dr_water, three_water &
    & , water_reach_ice, ice_water, three_rock &
    & , radius_ice_water, three_ice, frez_start, ice_rockice, ice_rock &
    & , material_select, ice_reach_rock &
    & , mass_ratio_rock, mass_ratio_anhydrous_rock, mass_ratio_hydrous_rock &
    & , mass_ratio_water, mass_ratio_ice &
    & , start_flux_alt, start_water_alt &
    & , mass_ratio_dehydration_rock, dehyd, mass_ratio_vapor, i_irw
  integer :: mn(0:irmax), irw, ir, iw, ir_old, iw_old, ii, ii_old, i_irw 
  real(rkind) :: dr_rock, radius_rock &
    & , radius_water, radius_ice, drock, dr_water, radius_ice_water &
    & , mass_ratio_rock(0:irmax), mass_ratio_anhydrous_rock(0:irmax) &
    & , mass_ratio_hydrous_rock(0:irmax), mass_ratio_water(0:irmax) &
    & , mass_ratio_ice(0:irmax), drock_wri, mass_ratio_dehydration_rock(0:irmax) &
    & , mass_ratio_vapor(0:irmax)
  logical :: melt(0:irmax) &
   & , water_reach_ice, frez_start, ice_reach_rock&
   & , start_flux_alt, start_water_alt, dehyd(0:irmax)
contains
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine material_init
    use global_constants, only : volume_ratio_rock, density &
      & , imr, mass_ratio_rock_H2O_init, itmax, Ordinary_Chondrite
    use global_parameters, only : radius
    real(rkind) :: mass_ratio_rock_mn6
    water_reach_ice=.FALSE.
    frez_start=.FALSE.
    ice_reach_rock=.FALSE.
    start_flux_alt=.FALSE.
    start_water_alt=.FALSE.
    if(.NOT.Ordinary_Chondrite) then
      mn(0:irmax) = 6 
      mass_ratio_rock_mn6 = mass_ratio_rock_H2O_init(imr)
      mass_ratio_rock(0:irmax) = mass_ratio_rock_mn6
      mass_ratio_anhydrous_rock(0:irmax) = mass_ratio_rock_mn6
      mass_ratio_ice(0:irmax) = 1._rkind - mass_ratio_rock_mn6
    elseif(Ordinary_Chondrite) then
      mn(0:irmax) = 4 
      mass_ratio_rock(0:irmax) = 1._rkind
      mass_ratio_anhydrous_rock(0:irmax) = 1._rkind
      mass_ratio_ice(0:irmax) = 0._rkind 
    endif
    mass_ratio_hydrous_rock(0:irmax) = 0._rkind
    mass_ratio_dehydration_rock(0:irmax) = 0._rkind
    mass_ratio_water(0:irmax) = 0._rkind
    mass_ratio_vapor(0:irmax) = 0._rkind
    melt(0:irmax) = .FALSE.
    dehyd(0:irmax) = .FALSE.
    ii = 0
    iw = 0
    ir = 0
    ir_old = 0
    iw_old = 0
    ii_old = 0
    radius_rock = 0._rkind
    radius_water = 0._rkind
    radius_ice = radius(0)
    drock = 0._rkind
  end subroutine material_init
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine material_select(i)
  integer,intent(in) :: i
    select case (mn(i))
      case (7) 
        call ice_water_rock(i)
      case (3) 
        call ice_water(i)
      case (8) 
        call three_water(i)
      case (5) 
        call rock_water(i)
      case (9)
        call three_ice(i)
      case (11)
        call water_rockice(i)
      case (12)
        call three_rock(i)
      case (13)
        call ice_rock(i)
      case (17)
        call rock_vapor_water(i)
    end select
  end subroutine material_select
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine material_check(i)
    use global_constants, only : rkind, volume_ratio_rock, time, year
    use global_parameters, only : temp, radius, halfdr
    use ice_parameters, only : temperature_melting_I
    integer,intent(in) :: i
    real(rkind) :: tmelt
    tmelt = temperature_melting_I
    if(temp(i) > tmelt) then
      melt(i) = .TRUE.
      call melt_check(i)
    elseif(temp(i) < tmelt) then
      melt(i) = .FALSE.
      if(mn(i)==7 .OR. mn(i) == 15) call frez_check(i)
    elseif(temp(i) == tmelt) then 
      if(melt(i)) call melt_check(i)
      if(.NOT. melt(i)) then
        if(mn(i)==7 .OR. mn(i) == 15) call frez_check(i)
      endif
    else
      call log_error(i,'material_check')
      return
    endif
    return
  end subroutine material_check
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine ice_water_rock(i)
    use global_constants, only : volume_ratio_rock, density
    use heat_parameters, only : thermal_conductivity, specific_heat & 
      & , thermometric_diffusivity
    integer,intent(in) :: i
    real(rkind) :: volume_ratio_ice, volume_ratio_water
    density(7) = 1._rkind/(mass_ratio_rock(i)/density(4) &
      & +mass_ratio_ice(i)/density(2)+mass_ratio_water(i)/density(1))
    volume_ratio_ice = mass_ratio_ice(i)*density(7)/density(2)
    volume_ratio_water = mass_ratio_water(i)*density(7)/density(1)
    thermal_conductivity(7) = thermal_conductivity(4)*volume_ratio_rock(7) &
      & + thermal_conductivity(2)*volume_ratio_ice &
      & + thermal_conductivity(1)*volume_ratio_water
    specific_heat(7) = specific_heat(4)*mass_ratio_rock(i) &
      & + specific_heat(2)*mass_ratio_ice(i)&
      & + specific_heat(1)*mass_ratio_water(i)
    thermometric_diffusivity(7) = thermal_conductivity(7) & 
      & /(density(7)*specific_heat(7))
  end subroutine ice_water_rock
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine rock_water(i)
    use global_constants, only : volume_ratio_rock, density
    use heat_parameters, only : thermal_conductivity, specific_heat & 
      & , thermometric_diffusivity
    integer,intent(in) :: i
    density(5) = volume_ratio_rock(5)*density(4) &
      & +(1._rkind-volume_ratio_rock(5))*density(1)
    mass_ratio_rock(i) = volume_ratio_rock(5)*density(4)/density(5)
    mass_ratio_water(i) = 1._rkind - mass_ratio_rock(i)
    thermal_conductivity(5) = thermal_conductivity(4)*volume_ratio_rock(5) &
      & + thermal_conductivity(1)*(1._rkind-volume_ratio_rock(5))
    specific_heat(5) = specific_heat(4)*mass_ratio_rock(i) &
      & + specific_heat(1)*(1._rkind-mass_ratio_rock(i))
    thermometric_diffusivity(5) = thermal_conductivity(5) & 
      & /(density(5)*specific_heat(5))
  end subroutine rock_water
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine three_water(i)
    use global_constants, only : volume_ratio_rock, density
    use global_parameters, only : radius, halfdr, dr
    use heat_parameters, only : thermal_conductivity, specific_heat & 
      & , thermometric_diffusivity
    integer,intent(in) :: i
    real(rkind) :: volume_ratio_water, mass_ratio_water8
    call ice_water_rock(i)
    volume_ratio_water = (radius_water-(radius(iw)+halfdr))/dr
    volume_ratio_rock(8) = (1._rkind-volume_ratio_water)*volume_ratio_rock(7)
    density(8) = density(1)*volume_ratio_water &
      & + density(7)*(1._rkind-volume_ratio_water)
    mass_ratio_water8 = volume_ratio_water*density(1)/density(8)
    thermal_conductivity(8) = thermal_conductivity(1)*volume_ratio_water &
      & + thermal_conductivity(7)*(1._rkind-volume_ratio_water)
    specific_heat(8) = specific_heat(1)*mass_ratio_water8 &
      & + specific_heat(7)*(1._rkind-mass_ratio_water8)
    thermometric_diffusivity(8) = thermal_conductivity(8) & 
      & /(density(8)*specific_heat(8))
  end subroutine three_water
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine ice_water(i)
    use global_constants, only : volume_ratio_rock, density
    use global_parameters, only : radius, halfdr, dr
    use heat_parameters, only : thermal_conductivity, specific_heat & 
      & , thermometric_diffusivity
    integer,intent(in) :: i
    real(rkind) :: volume_ratio_water
    mass_ratio_rock(i) = 0._rkind
    volume_ratio_water = (radius_water-(radius(iw)+halfdr))/dr
    density(3) = density(1)*volume_ratio_water &
      & + density(2)*(1._rkind-volume_ratio_water)
    mass_ratio_water(i) = volume_ratio_water*density(1)/density(3)
    mass_ratio_ice(i) = (1._rkind-volume_ratio_water)*density(2)/density(3)
    thermal_conductivity(3) = thermal_conductivity(1)*volume_ratio_water &
      & + thermal_conductivity(2)*(1._rkind-volume_ratio_water)
    specific_heat(3) = specific_heat(1)*mass_ratio_water(i) &
      & + specific_heat(2)*mass_ratio_ice(i)
    thermometric_diffusivity(3) = thermal_conductivity(3) & 
      & /(density(3)*specific_heat(3))
  end subroutine ice_water
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine three_ice(i)
    use global_constants, only : volume_ratio_rock, density
    use global_parameters, only : radius, halfdr, dr
    use heat_parameters, only : thermal_conductivity, specific_heat & 
      & , thermometric_diffusivity
    integer,intent(in) :: i
    real(rkind) :: volume_ratio_ice, volume_ratio_water, volume_ratio_H2O
    volume_ratio_ice = (radius_ice_water-radius_water)/dr
    volume_ratio_water = (radius_water-(radius(iw)+halfdr))/dr
    volume_ratio_H2O = volume_ratio_ice + volume_ratio_water
    volume_ratio_rock(9) = (1._rkind-volume_ratio_H2O)*volume_ratio_rock(6)
    density(9) = density(2)*volume_ratio_ice + density(1)*volume_ratio_water &
      & + density(6)*(1._rkind-volume_ratio_H2O)
    mass_ratio_ice(i) = (volume_ratio_ice+(1-volume_ratio_H2O)&
      & *(1._rkind-volume_ratio_rock(6)))*density(2)/density(9)
    mass_ratio_water(i) = volume_ratio_water*density(1)/density(9)
    mass_ratio_rock(i) = volume_ratio_rock(9)*density(4)/density(9)
    thermal_conductivity(9) = thermal_conductivity(2)*volume_ratio_ice &
      & + thermal_conductivity(1)*volume_ratio_water &
      & + thermal_conductivity(4)*volume_ratio_rock(9)
    specific_heat(9) = specific_heat(2)*mass_ratio_ice(i) &
      & + specific_heat(1)*mass_ratio_water(i) &
      & + specific_heat(4)*mass_ratio_rock(i)
   thermometric_diffusivity(9) = thermal_conductivity(9) & 
      & /(density(9)*specific_heat(9))
  end subroutine three_ice
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine water_rockice(i)
    use global_constants, only : volume_ratio_rock, density
    use global_parameters, only : radius, halfdr, dr
    use heat_parameters, only : thermal_conductivity, specific_heat & 
      & , thermometric_diffusivity
    integer,intent(in) :: i
    real(rkind) :: volume_ratio_water, mass_ratio_water11
    volume_ratio_water = (radius_water-(radius(i)-halfdr))/dr
    volume_ratio_rock(11) = (1._rkind-volume_ratio_water)*volume_ratio_rock(6)
    density(11) = density(1)*volume_ratio_water + density(6)*(1._rkind-volume_ratio_water)
    mass_ratio_water11 = volume_ratio_water*density(1)/density(11)
    thermal_conductivity(11) = thermal_conductivity(1)*volume_ratio_water &
      & + thermal_conductivity(6)*(1._rkind-volume_ratio_water)
    specific_heat(11) = specific_heat(1)*mass_ratio_water11 &
      & + specific_heat(6)*(1._rkind-mass_ratio_water11)
   thermometric_diffusivity(11) = thermal_conductivity(11) & 
      & /(density(11)*specific_heat(11))
  end subroutine water_rockice
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine ice_rockice(i)
    use global_constants, only : volume_ratio_rock, density
    use global_parameters, only : radius, halfdr, dr
    use heat_parameters, only : thermal_conductivity, specific_heat & 
      & , thermometric_diffusivity
    integer,intent(in) :: i
    real(rkind) :: volume_ratio_ice
    volume_ratio_ice = (radius_ice_water-(radius(ii)-halfdr))/dr
    volume_ratio_rock(10) = (1._rkind-volume_ratio_ice)*volume_ratio_rock(6)
    density(10) = density(2)*volume_ratio_ice + density(6)*(1._rkind-volume_ratio_ice)
    mass_ratio_water(i) = 0._rkind
    mass_ratio_ice(i) = (volume_ratio_ice+(1-volume_ratio_ice)&
      & *(1._rkind-volume_ratio_rock(6)))*density(2)/density(10)
    mass_ratio_rock(i) = volume_ratio_rock(10)*density(4)/density(10)
    thermal_conductivity(10) = thermal_conductivity(2)*volume_ratio_ice &
      & + thermal_conductivity(6)*(1._rkind-volume_ratio_ice)
    specific_heat(10) = specific_heat(2)*mass_ratio_ice(i) &
      & + specific_heat(6)*(1._rkind-mass_ratio_ice(i))
   thermometric_diffusivity(10) = thermal_conductivity(10) & 
      & /(density(10)*specific_heat(10))
  end subroutine ice_rockice
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine three_rock(i)
    use global_constants, only : volume_ratio_rock, density
    use global_parameters, only : radius, halfdr, dr
    use heat_parameters, only : thermal_conductivity, specific_heat & 
      & , thermometric_diffusivity
    integer,intent(in) :: i
    real(rkind) :: volume_ratio_ice, volume_ratio_water
    volume_ratio_ice = ((radius(ii)+halfdr)-radius_water)/dr
    volume_ratio_water = (radius_water-radius_rock)/dr
    volume_ratio_rock(12) = (radius_rock-(radius(ii-1)+halfdr))/dr
    density(12) = density(4)*volume_ratio_rock(12) &
      & + density(1)*volume_ratio_water + density(2)*volume_ratio_ice
    mass_ratio_ice(i) = volume_ratio_ice*density(2)/density(12)
    mass_ratio_water(i) = volume_ratio_water*density(1)/density(12)
    mass_ratio_rock(i) = volume_ratio_rock(12)*density(4)/density(12)
    thermal_conductivity(12) = thermal_conductivity(4)*volume_ratio_rock(12) &
      & + thermal_conductivity(2)*volume_ratio_ice &
      & + thermal_conductivity(1)*volume_ratio_water 
    specific_heat(12) = specific_heat(4)*mass_ratio_rock(i) &
      & + specific_heat(2)*mass_ratio_ice(i) &
      & + specific_heat(1)*mass_ratio_water(i) 
    thermometric_diffusivity(12) = thermal_conductivity(12) & 
      & /(density(12)*specific_heat(12))
  end subroutine three_rock
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine ice_rock(i)
    use global_constants, only : volume_ratio_rock, density
    use global_parameters, only : radius, halfdr, dr
    use heat_parameters, only : thermal_conductivity, specific_heat & 
      & , thermometric_diffusivity
    integer,intent(in) :: i
    real(rkind) :: volume_ratio_ice
    volume_ratio_ice = ((radius(ii)+halfdr)-radius_rock)/dr
    volume_ratio_rock(13) = (radius_rock-(radius(ii-1)+halfdr))/dr
    density(13) = density(4)*volume_ratio_rock(13) + density(2)*volume_ratio_ice
    mass_ratio_water(i) = 0._rkind
    mass_ratio_ice(i) = volume_ratio_ice*density(2)/density(13)
    mass_ratio_rock(i) = volume_ratio_rock(13)*density(4)/density(13)
    thermal_conductivity(13) = thermal_conductivity(4)*volume_ratio_rock(13) &
      & + thermal_conductivity(2)*volume_ratio_ice 
    specific_heat(13) = specific_heat(4)*mass_ratio_rock(i) &
      & + specific_heat(2)*mass_ratio_ice(i)
    thermometric_diffusivity(13) = thermal_conductivity(13) & 
      & /(density(13)*specific_heat(13))
  end subroutine ice_rock
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine rock_vapor_water(i)
    use global_constants, only : volume_ratio_rock, density
    use heat_parameters, only : thermal_conductivity, specific_heat & 
      & , thermometric_diffusivity
    integer,intent(in) :: i
    real(rkind) :: volume_ratio_vapor, volume_ratio_water
    density(17) = 1._rkind/(mass_ratio_rock(i)/density(4) &
      & +mass_ratio_vapor(i)/density(0)+mass_ratio_water(i)/density(1))
    volume_ratio_vapor = mass_ratio_vapor(i)*density(17)/density(0)
    volume_ratio_water = mass_ratio_water(i)*density(17)/density(1)
    thermal_conductivity(17) = thermal_conductivity(4)*volume_ratio_rock(17) &
      & + thermal_conductivity(0)*volume_ratio_vapor &
      & + thermal_conductivity(1)*volume_ratio_water
    specific_heat(17) = specific_heat(4)*mass_ratio_rock(i) &
      & + specific_heat(0)*mass_ratio_vapor(i)&
      & + specific_heat(1)*mass_ratio_water(i)
    thermometric_diffusivity(17) = thermal_conductivity(17) & 
      & /(density(17)*specific_heat(17))
  end subroutine rock_vapor_water
!--------1---------2---------3---------4---------5---------6---------7---------8
end module material_parameters