!--------1---------2---------3---------4---------5---------6---------7---------8
module global_constants
!--------1---------2---------3---------4---------5---------6---------7---------8
  implicit none
  public
  integer,parameter :: rkind = kind(0.d0)
  integer,parameter :: qkind = kind(0.q0)
  integer,parameter :: itmax = 50000000  
  integer,parameter :: cat = 30
  integer,parameter :: istep_max = 100 
  character(len=cat),save :: ice_name = 'IceOnly'
  real(rkind),parameter :: emissivity = 1._rkind
  real(rkind),parameter :: absorption_coefficient = 1._rkind
!--------1---------2---------3---------4---------5---------6---------7---------8
  real(rkind),parameter :: pi = 3.14159265358979323846264338327950_rkind
  real(rkind),parameter :: stefan_boltzman_constant = 5.67e-8_rkind
  real(rkind),parameter :: luminosity_star = 3.86e26_rkind
  real(rkind),parameter :: mass_Hatom = 1.673e-27_rkind
  real(rkind),parameter :: MeV_J = 1.6e-13_rkind
  real(rkind),parameter :: gravitational_constant = 6.672e-11_rkind
  real(rkind),parameter :: micrometer = 1.e-6_rkind
  real(rkind),parameter :: kilometer = 1e3_rkind
  real(rkind),parameter :: au = 1.496e8*kilometer
  real(rkind),parameter :: temperature_Kelvin = 273.16_rkind
  real(rkind),parameter :: thickness_epsilon = 1._rkind
  real(rkind),parameter :: minute = 60._rkind
  real(rkind),parameter :: hour = 60._rkind*minute
  real(rkind),parameter :: day = 24._rkind*hour
  real(rkind),parameter :: year = 365._rkind*day
  real(rkind),parameter :: Myr = 1.e6_rkind
  real(rkind),parameter :: time0 = 1._rkind
  real(rkind),parameter :: timeL = 1.e7_rkind*year
  real(rkind),parameter :: dt = (timeL - time0)/itmax
  real(rkind) :: time, density(0:20), volume_ratio_rock(0:20), time_aqueous
  real(rkind),parameter :: density_vapor = 1.e3_rkind
  real(rkind),parameter :: density_water = 1.e3_rkind
  real(rkind),parameter :: density_ice0 = 1.e3_rkind
  real(rkind),parameter :: density_ice = 0.93e3_rkind 
  real(rkind),parameter :: density_rock = 1.82e3_rkind 
  real(rkind),parameter :: density_rock_hydrous = 1.1_rkind*density_rock*density_water/(density_water+0.1_rkind*density_rock)
  real(rkind),parameter :: volume_ratio_rock_pure = 1._rkind
  real(rkind),parameter :: volume_ratio_rock_non = 0._rkind
  real(rkind),parameter :: mass_ratio_rock_H2O_solar = 0.625_rkind 
  real(rkind),parameter :: temperature_melting_rock = 1773.16_rkind
  real(rkind),parameter :: temperature_aqueous_alteration = 293.16_rkind
  real(rkind),parameter :: cool_temperature = 400._rkind + temperature_Kelvin
  real(rkind),parameter :: cool_temperature_TOP = cool_temperature + 0._rkind
  real(rkind),parameter :: cool_temperature_BOTTOM = cool_temperature - 1._rkind
  real(rkind),parameter :: mass_ratio_rock_new = 0.82_rkind 
  real(rkind),parameter :: specific_heat_rockA = 564._rkind 
  real(rkind),parameter :: specific_heat_rockB = 0._rkind
  real(rkind),parameter :: specific_heat_rockC = 0._rkind
  real(rkind),parameter :: thermometric_conductivity_rockA = 3.2e-7_rkind 
  real(rkind),parameter :: thermometric_conductivity_rockB = 0._rkind
  real(rkind),parameter :: specific_heat_ice0 = 1925._rkind
  real(rkind),parameter :: specific_heat_iceA = 1380._rkind 
  real(rkind),parameter :: thermometric_conductivity_iceA = 3.69_rkind/specific_heat_iceA/density_ice 
  real(rkind),parameter :: thermometric_conductivity_ice0 = 2.2_rkind/specific_heat_ice0/density_ice0
  real(rkind),parameter :: specific_heat_waterA = 4200._rkind
  real(rkind),parameter :: thermometric_conductivity_waterA = 0.56_rkind/specific_heat_waterA/density_water
  logical :: GO_NEXT = .FALSE., W_ALTERATION, D_REACTION, WV_PRESSURE &
  & , Rocky_Core, Ordinary_Chondrite, Aqueous_Start = .FALSE.
  integer :: imr
  real(rkind) :: mass_ratio_rock_H2O_init(1:9)
contains
  subroutine set_mass_ratio_rock_H2O_init
    mass_ratio_rock_H2O_init(1) = mass_ratio_rock_H2O_solar
    do imr = 2, 9
      mass_ratio_rock_H2O_init(imr) = real(imr)
    enddo
  end subroutine set_mass_ratio_rock_H2O_init
  subroutine set_density
    density(0) = density_vapor
    density(1) = density_water
    density(2) = density_ice
    density(4) = density_rock
    density(6) = 1._rkind/(mass_ratio_rock_H2O_init(imr)/density_rock &
      & + (1._rkind-mass_ratio_rock_H2O_init(imr))/density_ice)
    density(15) = 1._rkind/(mass_ratio_rock_new/density_rock &
      & + (1._rkind-mass_ratio_rock_new)/density_water)
    density(16) = 1._rkind/(mass_ratio_rock_new/density_rock &
      & + (1._rkind-mass_ratio_rock_new)/density_vapor)
  end subroutine set_density
  subroutine set_volume_ratio_rock
    real(rkind) :: volume_ratio_rock_H2O
    volume_ratio_rock_H2O = mass_ratio_rock_H2O_init(imr)*density(6)/density_rock
    volume_ratio_rock(1) = volume_ratio_rock_non
    volume_ratio_rock(2) = volume_ratio_rock_non
    volume_ratio_rock(3) = volume_ratio_rock_non
    volume_ratio_rock(4) = volume_ratio_rock_pure
    volume_ratio_rock(5) = volume_ratio_rock_H2O
    volume_ratio_rock(6) = volume_ratio_rock_H2O
    volume_ratio_rock(7) = volume_ratio_rock_H2O
    volume_ratio_rock(15) = mass_ratio_rock_new*density(15)/density_rock
    volume_ratio_rock(16) = mass_ratio_rock_new*density(16)/density_rock
  end subroutine set_volume_ratio_rock
end module global_constants
!--------1---------2---------3---------4---------5---------6---------7---------8
module global_parameters 
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, pi, gravitational_constant, density, imr
  implicit none
  private
  public ::  set_temperature_initial, temperature_initial, iau, iaumax & 
    & , set_radius_body, radius_body, irb, irbmax &
    & , set_radius, radius, dr, set_temperature, temp, irmax &
    & , pressure_ice & 
    & , surface_area, surface_arear, halfdr, lat_heat, rad_heat &
    & , fw, fname_center, fname_all, fname_E, fname_Eall, fname_Enew, fname_Ecum &
    & , pressure_vapor, pressure_lithostatic &
    & , fname_radius, temp_max, cooling_rate, cooled_time &
    & , igmax, concentration, diffusion_length, fname_max
  integer,parameter :: irbmax = 200
  integer,parameter :: iaumax = 2 
  integer,parameter :: irmax = 400
  integer,parameter :: igmax = 100
  integer :: irb, iau, fw
  real(rkind) :: radius_body(1:irbmax), temperature_initial(1:iaumax), halfdr &
    & , radius(0:irmax), temp(0:irmax), dr, acceleration_gravity(1:irbmax) &
    & , lat_heat(0:irmax), rad_heat(0:irmax), pressure_lithostatic(0:irmax) &
    & , temp_max(0:irmax), cooling_rate(0:irmax), cooled_time(0:irmax) &
    & , concentration(0:irmax,0:igmax), diffusion_length(0:irmax)
  character*16 :: fname_center, fname_all, fname_E, fname_Eall, fname_Enew, fname_Ecum, fname_radius, fname_max
contains
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine set_radius_body
    use global_constants, only : kilometer
    do irb = 1, irbmax
      radius_body(irb) = kilometer*real(irb)
    enddo
    return
  end subroutine set_radius_body
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine set_temperature_initial
    use global_constants, only : au, luminosity_star &
      & , stefan_boltzman_constant, absorption_coefficient, emissivity
     temperature_initial(1) = 70._rkind
     temperature_initial(2) = 150._rkind
    return
  end subroutine set_temperature_initial
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine set_radius
    integer :: ir
    dr = radius_body(irb)/real(irmax)
    halfdr = dr*0.5_rkind
    do ir = 0, irmax
      radius(ir) = dr*real(ir)
    enddo
  return
  end subroutine set_radius
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine set_temperature
    integer :: ir
    do ir = 0, irmax
      temp(ir) = temperature_initial(iau)
      temp_max(ir) = temperature_initial(iau)
      cooling_rate(ir) = 0._rkind
      cooled_time(ir) = 0._rkind
    enddo
    return
  end subroutine set_temperature
!--------1---------2---------3---------4---------5---------6---------7---------8
  real(rkind) function pressure_ice(thick)
    real(rkind),intent(in) :: thick
    real(rkind) :: coef_pressure
    coef_pressure = 2._rkind/3._rkind*pi*density(6)*density(6) &
      & *gravitational_constant 
    pressure_ice = coef_pressure*(2._rkind*radius_body(irb)-thick)*thick
  end function pressure_ice
  real(rkind) function surface_area(i)
    implicit none
    integer,intent(in) :: i
    surface_area = 4._rkind*pi*radius(i)*radius(i)
    return
  end function surface_area
  real(rkind) function surface_arear(r)
    implicit none
    real(rkind),intent(in) :: r
    surface_arear = 4._rkind*pi*r*r
    return
  end function surface_arear
!--------1---------2---------3---------4---------5---------6---------7---------8
  real(rkind) function pressure_vapor(temp)
    real(rkind),intent(in) :: temp
    pressure_vapor = 100._rkind*exp(-6096.9385_rkind/temp+16.635794 &
      & -2.711193e-2_rkind*temp+1.673952e-5_rkind*temp*temp &
      & +2.433502_rkind*log(temp))
  end function pressure_vapor
!--------1---------2---------3---------4---------5---------6---------7---------8
end module global_parameters