!--------1---------2---------3---------4---------5---------6---------7---------8
module aqueous 
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind
  use global_parameters, only : irmax
  implicit none
  private
  public :: aqueous_check, aqueous_alteration, aqueous_alteration_below &
    & , alteration_heat_bak, flux_alt
  real(rkind),parameter :: alt_ratio_water_anhydrous = 0.1_rkind
  real(rkind),parameter :: alt_ratio_hydrous_anhydrous = 1.1_rkind
  real(rkind),parameter :: alteration_heat = 2.77e5_rkind
  real(rkind) :: d_hydrous(0:irmax), flux_alt, drw
contains
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine aqueous_check
    use global_constants, only : temperature_aqueous_alteration
    use global_parameters, only : temp
    use material_parameters, only : mn, mass_ratio_anhydrous_rock
    integer :: i
    do i = 0, irmax
      if(mass_ratio_anhydrous_rock(i) <= 0._rkind) cycle
      if(temp(i) >= temperature_aqueous_alteration) then
        if(mn(i) == 8 .OR. mn(i) == 7 .OR. mn(i) == 15) then
          call aqueous_alteration(i)
        endif
      endif
    enddo
  end subroutine aqueous_check
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine aqueous_alteration(i)
    use global_constants, only : time, year, Myr, dt, temperature_aqueous_alteration
    use global_parameters, only : pressure_lithostatic
    use heat_radiogenic, only : iAl, ftime
    use material_parameters, only : mass_ratio_water, mn &
      & , mass_ratio_rock, mass_ratio_anhydrous_rock, mass_ratio_hydrous_rock, mass_ratio_ice
    integer,intent(in) :: i
    real(rkind) :: need_water, d_anhydrous, d_water, time_Myr &
      & , reaction_rate, aqueous_rate 
    mass_ratio_hydrous_rock(i) = 1.1_rkind*mass_ratio_anhydrous_rock(i)                  
    d_hydrous(i) = mass_ratio_hydrous_rock(i)
    mass_ratio_water(i) = 1._rkind - mass_ratio_hydrous_rock(i)
    mass_ratio_anhydrous_rock(i) = 0._rkind 
    mass_ratio_ice(i) = 0._rkind
    mass_ratio_rock(i) = mass_ratio_hydrous_rock(i) + mass_ratio_anhydrous_rock(i)
    time_Myr = (time/year+ftime(iAl))/Myr
  end subroutine aqueous_alteration
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine aqueous_alteration_below(i,dr_water)
    use global_constants, only : density, time, year
    use material_parameters, only : mass_ratio_ice, mass_ratio_anhydrous_rock &
      & , mass_ratio_hydrous_rock, drock, radius_water, mn
    integer,intent(in) :: i
    real(rkind),intent(in) :: dr_water
    real(rkind) :: need_water, d_anhydrous, d_water, mrw, density_tmp
    if(mass_ratio_anhydrous_rock(i) > 0._rkind) then
      need_water = alt_ratio_water_anhydrous*mass_ratio_anhydrous_rock(i)
      mrw = mass_ratio_ice(i)
      if(mrw >= need_water) then
        d_water = need_water
        d_anhydrous = mass_ratio_anhydrous_rock(i)
      elseif(mrw < need_water) then
        d_water = mrw
        d_anhydrous = d_water/alt_ratio_water_anhydrous
      endif
      d_hydrous(i) = d_hydrous(i) + alt_ratio_hydrous_anhydrous*d_anhydrous
      d_anhydrous = mass_ratio_anhydrous_rock(i) - d_anhydrous
      density_tmp = 1._rkind/(d_hydrous(i)/density(4)+d_anhydrous/density(4) &
        & +(mrw-d_water)/density(1))
      drock = drock + (d_hydrous(i)/density(4)+d_anhydrous/density(4)) &
        & *density_tmp*dr_water*radius_water*radius_water
      drw = dr_water
    elseif(mass_ratio_anhydrous_rock(i) <= 0._rkind) then
      drock = drock + mass_ratio_hydrous_rock(i)/density(4) &
        & *density(6)*dr_water*radius_water*radius_water
      return
    endif
  end subroutine aqueous_alteration_below
!--------1---------2---------3---------4---------5---------6---------7---------8
  subroutine alteration_heat_bak(it)
    use global_constants, only : imr, density, dt, time, year, pi, Aqueous_Start, time_aqueous &
      & , mass_ratio_rock_H2O_init, Myr, density_rock_hydrous, mass_ratio_rock_new, mass_ratio_rock_H2O_solar
    use global_parameters, only : surface_area, dr, temp, radius, halfdr
    use heat_parameters, only : latent_heat_init, specific_heat
    use material_parameters, only : mn, mass_ratio_ice, mass_ratio_water &
      & , radius_water, start_flux_alt, start_water_alt, mass_ratio_rock, material_check &
      & , mass_ratio_anhydrous_rock, mass_ratio_hydrous_rock, mass_ratio_dehydration_rock
    use energy, only : lat_energy, alt_energy, alt_lat_energy, temp_old
    use energy_dump, only : it_start_water_alt, it_start_flux_alt
  use heat_radiogenic, only : iAl, ftime
    integer,intent(in) :: it
    integer :: i
    real(rkind) :: dm, derg, dm_new, time_Myr, dtemp, Tspecific_heat, tmp, Tdensity
    time_Myr = (time/year+ftime(iAl))/Myr
    do i = 0, irmax
      if(d_hydrous(i) == 0._rkind) cycle
      if(mn(i) == 7 .OR. mn(i) == 8 .OR. mn(i) == 15) then
          if(.NOT.Aqueous_Start) then
            Aqueous_Start=.TRUE.
            time_aqueous = time_Myr
          endif
        dm = alteration_heat*d_hydrous(i)/latent_heat_init
        if(i/=0) then
          derg = latent_heat_init*dm*Tdensity(i)*surface_area(i)*dr/dt
        elseif(i==0) then
          derg = latent_heat_init*dm*Tdensity(i)*4._rkind/3._rkind*pi*dr*dr*dr/dt
        endif
        if(.NOT.start_water_alt) then
          start_water_alt=.TRUE.
          it_start_water_alt=it
        endif
          dtemp = (alteration_heat*d_hydrous(i) - (1._rkind-mass_ratio_rock_H2O_solar)*latent_heat_init)/Tspecific_heat(i)
          call material_check(i)
          dm_new = dm + mass_ratio_ice(i)
          temp(i) = temp(i) + dtemp
          alt_lat_energy(i) = derg*dm_new/dm
          alt_energy(i) = derg*dm_new/dm
          mass_ratio_ice(i) = 0._rkind
          mass_ratio_water(i) = 1._rkind - mass_ratio_rock(i)
          radius_water = radius(i) + halfdr
      elseif(mn(i) == 11 .OR. mn(i) == 1) then
        derg = alteration_heat*d_hydrous(i)*Tdensity(i)*surface_area(i)*drw/dt
        if(.NOT.start_flux_alt) then
          start_flux_alt=.TRUE.
          it_start_flux_alt=it
        endif
        flux_alt = flux_alt + derg
        alt_energy(i) = derg
      elseif(mn(i) == 15) then
        write(1,'(f11.7,2i4,3e12.5,a)') time_Myr, i, mn(i)&
          & , mass_ratio_hydrous_rock(i), mass_ratio_anhydrous_rock(i), mass_ratio_water(i), 'alt15' 
        derg = alteration_heat*d_hydrous(i)
        dtemp = (derg - (1._rkind-mass_ratio_rock_H2O_solar)*latent_heat_init)/Tspecific_heat(i)
        temp(i) = temp(i) + dtemp
        if(i/=0) then
          alt_energy(i) = derg*Tdensity(i)*surface_area(i)*dr/dt
        elseif(i==0) then
          alt_energy(i) = derg*Tdensity(i)*4._rkind/3._rkind*pi*dr*dr*dr/dt
        endif
      else
        write(999,*) 'Sorry this program is end here'
        write(999,*) i, mn(i), d_hydrous(i)
        call log_error(i,'alteration_heat_bak mn')
        stop
      endif
    enddo
    d_hydrous(0:irmax) = 0._rkind   
  end subroutine alteration_heat_bak
end module aqueous