!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine pressure_check
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, gravitational_constant, pi, dt &
    & , density, kilometer, WV_PRESSURE, time, year, Myr
  use heat_radiogenic, only : iAl, ftime
  use global_parameters, only : irmax, radius, temp, dr, surface_area &
    & , pressure_vapor, pressure_lithostatic
  use heat_parameters, only : specific_heat, latent_heat_vapor
  use material_parameters, only : radius_water, radius_rock, ir, iw, i_irw &
    & , mass_ratio_vapor, mass_ratio_water, material_select, mn &
    & , mass_ratio_rock, frez_start, radius_ice_water
  use energy, only : lat_energy, temp_old, rad_energy
  implicit none
  integer :: i, i_mantle, i_core
  real(rkind) :: coef, coef_crust, coef_mantle, radius_mantle, radius_core &
    & , m, dm, coef_dm, latent_heat, temp_new, time_Myr, g(0:irmax)
  coef = 2._rkind/3._rkind*gravitational_constant*pi
  coef_mantle = radius_rock**3._rkind*(density(15)-density(1)) &
    & *density(1)*2._rkind*coef
  if(frez_start) then
    i_mantle = i_irw
    i_core = ir
    radius_mantle = radius_ice_water
    radius_core = radius_rock
  elseif(.NOT.frez_start) then 
    i_mantle = iw
    i_core = ir
    radius_mantle = radius(iw)
    radius_core = radius(ir)
  endif
  do i = i_mantle+1, irmax-1
    pressure_lithostatic(i) = coef*density(6)*density(6) &
      & *(radius(irmax)+radius(i))*(radius(irmax)-radius(i)) 
    g(i) = coef*2._rkind*radius(i)*density(6)
  enddo
  do i = ir+1, i_mantle
    pressure_lithostatic(i) = pressure_lithostatic(i_mantle+1) + &
      & coef*density(1)*density(1)*(radius_mantle+radius(i))*(radius_mantle-radius(i)) &
      & + coef_mantle*(1._rkind/radius(i)-1._rkind/radius_mantle)
    g(i) = coef*2._rkind*(radius(i)*density(1) + radius_rock**3._rkind*(density(15)-density(1))/radius(i)/radius(i))
  enddo
  do i = 0, i_core
    pressure_lithostatic(i) = pressure_lithostatic(ir+1) + &
      & coef*density(15)*density(15)*(radius_core+radius(i))*(radius_core-radius(i)) 
    g(i) = coef*2._rkind*radius(i)*density(15)
  enddo
  time_Myr = (time/year+ftime(iAl))/Myr
  do i = 0, i_core
   if(WV_PRESSURE) then
   if(temp(i) >= temp_old(i)) then
     if(mass_ratio_water(i) /= 0._rkind .AND. pressure_vapor(temp(i)) >= pressure_lithostatic(i)) then
       m = mass_ratio_water(i)
       call material_select(i)
       latent_heat = latent_heat_vapor
       coef_dm = specific_heat(mn(i))/latent_heat
         dm = (temp(i)-temp_old(i))*coef_dm
         lat_energy(i) = lat_energy(i) + latent_heat*dm*density(mn(i))*surface_area(i)*dr/dt
         temp(i) = temp_old(i)
       m = m - dm
       mass_ratio_vapor(i) = mass_ratio_vapor(i) + dm 
       mn(i) = 17
       if(m <= 1e-7_rkind) then
         m = 0._rkind
         mn(i) = 16
       endif
       mass_ratio_water(i) = m
     endif
   elseif(temp(i) < temp_old(i)) then
     if(.NOT.frez_start) cycle
     if(mass_ratio_vapor(i) == 0._rkind) cycle
     if(mass_ratio_vapor(i) /= 0._rkind .AND. pressure_vapor(temp(i)) <= pressure_lithostatic(i)) then
         m = mass_ratio_vapor(i)
         call material_select(i)
         latent_heat = latent_heat_vapor
         coef_dm = specific_heat(mn(i))/latent_heat
       if(mass_ratio_water(i) /= 0._rkind) then
         dm = (temp_old(i)-temp(i))*coef_dm
         lat_energy(i) = lat_energy(i) + latent_heat*dm*density(mn(i))*surface_area(i)*dr/dt
         temp(i) = temp_old(i)
       elseif(mass_ratio_vapor(i) == 0._rkind) then
         call bisection_method(temp_new,pressure_lithostatic(i),i)
         dm = (temp_new-temp(i))*coef_dm
         lat_energy(i) = lat_energy(i) + latent_heat*dm*density(mn(i))*surface_area(i)*dr/dt
         temp(i) = temp_new
       endif
       m = m - dm
       mass_ratio_water(i) = mass_ratio_water(i) + dm 
       mn(i) = 17
       if(m <= 1e-6_rkind) then
         m = 0._rkind
         mn(i) = 15
       endif
       mass_ratio_vapor(i) = m
     endif
   endif
   endif
  enddo
end subroutine pressure_check
!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine pressure_check_nocore
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, gravitational_constant, pi, dt &
    & , density, kilometer, WV_PRESSURE, time, year, Myr, Rocky_Core
  use heat_radiogenic, only : iAl, ftime
  use global_parameters, only : irmax, radius, temp, dr, surface_area &
    & , pressure_vapor, pressure_lithostatic
  use heat_parameters, only : specific_heat, latent_heat_vapor
  use material_parameters, only : radius_water, radius_rock, ir, iw, i_irw, ii &
    & , mass_ratio_vapor, mass_ratio_water, material_select, mn &
    & , mass_ratio_rock, frez_start, radius_ice_water
  use energy, only : lat_energy, temp_old, rad_energy
  implicit none
  integer :: i, i_mantle, i_core
  real(rkind) :: coef, coef_crust, coef_mantle, radius_mantle, radius_core &
    & , m, dm, coef_dm, latent_heat, temp_new, time_Myr, g(0:irmax)
  coef = 2._rkind/3._rkind*gravitational_constant*pi
  coef_mantle = radius_rock**3._rkind*(density(15)-density(1)) &
    & *density(1)*2._rkind*coef
  radius_mantle = radius(irmax)
  do i = 0, irmax-1
    pressure_lithostatic(i) = coef*density(6)*density(6) &
      & *(radius(irmax)+radius(i))*(radius(irmax)-radius(i)) 
    g(i) = coef*2._rkind*radius(i)*density(6)
  enddo
  time_Myr = (time/year+ftime(iAl))/Myr
  do i = 0, irmax-1
   if(WV_PRESSURE) then
   if(temp(i) >= temp_old(i)) then
     if(mass_ratio_water(i) /= 0._rkind .AND. pressure_vapor(temp(i)) >= pressure_lithostatic(i)) then
       m = mass_ratio_water(i)
       call material_select(i)
       latent_heat = latent_heat_vapor
       coef_dm = specific_heat(mn(i))/latent_heat
         dm = (temp(i)-temp_old(i))*coef_dm
         lat_energy(i) = lat_energy(i) + latent_heat*dm*density(mn(i))*surface_area(i)*dr/dt
         temp(i) = temp_old(i)
       m = m - dm
       mass_ratio_vapor(i) = mass_ratio_vapor(i) + dm 
       mn(i) = 17
       if(m <= 1e-7_rkind) then
         m = 0._rkind
         mn(i) = 16
       endif
       mass_ratio_water(i) = m
     endif
   elseif(temp(i) < temp_old(i)) then
     if(mass_ratio_vapor(i) == 0._rkind) cycle
     if(mass_ratio_vapor(i) /= 0._rkind .AND. pressure_vapor(temp(i)) <= pressure_lithostatic(i)) then
         m = mass_ratio_vapor(i)
         call material_select(i)
         latent_heat = latent_heat_vapor
         coef_dm = specific_heat(mn(i))/latent_heat
       if(mass_ratio_water(i) /= 0._rkind) then
         dm = (temp_old(i)-temp(i))*coef_dm
         lat_energy(i) = lat_energy(i) + latent_heat*dm*density(mn(i))*surface_area(i)*dr/dt
         temp(i) = temp_old(i)
       elseif(mass_ratio_vapor(i) == 0._rkind) then
         call bisection_method(temp_new,pressure_lithostatic(i),i)
         dm = (temp_new-temp(i))*coef_dm
         lat_energy(i) = lat_energy(i) + latent_heat*dm*density(mn(i))*surface_area(i)*dr/dt
         temp(i) = temp_new
       endif
       m = m - dm
       mass_ratio_water(i) = mass_ratio_water(i) + dm 
       mn(i) = 17
       if(m <= 1e-6_rkind) then
         m = 0._rkind
         mn(i) = 15
       endif
       mass_ratio_vapor(i) = m
     endif
   endif
   endif
  enddo
end subroutine pressure_check_nocore
!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine bisection_method(temp_new,P_litho,i)
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, istep_max
  use global_parameters, only : temp, pressure_vapor
  use energy, only : temp_old
  implicit none
  real(rkind),intent(out) :: temp_new
  integer,intent(in) :: i
  real(rkind),intent(in) :: P_litho
  real(rkind) :: dmid, dstart, dend
  integer :: istep
  dstart = temp_old(i)
  dend = temp(i)
  do istep = 1, istep_max
    dmid = (dstart + dend)*0.5_rkind
    write(2,'(2i5,2f10.3,3e13.5)') i,istep,dstart,dend,pressure_vapor(dmid),P_litho,abs(pressure_vapor(dmid)-P_litho)
    if(abs(pressure_vapor(dmid) - P_litho) <= 1.e-1_rkind) exit
    if((pressure_vapor(dmid) - P_litho)*(pressure_vapor(dstart) - P_litho) <= 0._rkind ) then
      dend = dmid
    else
      dstart = dmid
    endif
    if(istep >= istep_max) then
      write(999,*) 'istep OVER istep_max'
      stop
    endif
  enddo
  temp_new = dmid
  return
end subroutine bisection_method