!--------1---------2---------3---------4---------5---------6---------7---------8
module diffusion
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : rkind, micrometer, temperature_Kelvin
  use global_parameters, only : irmax, igmax, concentration, diffusion_length
  implicit none
  private
  public :: diffusion_init, diffusion_check
  real(rkind),parameter :: gas_constant = 8.3144621_rkind
  real(rkind),parameter :: gas_constant_kcal = 1.987e-3_rkind
!--------1---------2---------3---------4---------5---------6---------7---------8
  real(rkind),parameter :: all_distance = 1._rkind*micrometer
  real(rkind),parameter :: dg = all_distance/real(igmax)
  real(rkind),parameter :: concentration_initial = 100._rkind
  real(rkind),parameter :: concentration_boundary = 0._rkind
  real(rkind),parameter :: temperature_boundary = 500._rkind+temperature_Kelvin
  real(rkind),parameter :: diff_eps = 1.e-7_rkind
contains
!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine diffusion_init
  diffusion_length(0:irmax) = 0._rkind
end subroutine diffusion_init
!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine diffusion_check(time,dt,temp,temp_old) 
!--------1---------2---------3---------4---------5---------6---------7---------8
  real(rkind),intent(in) :: time, dt, temp(0:irmax), temp_old(0:irmax)
  integer :: i
  do i = 0, irmax
        call calc_diffusion_length(time,dt,temp(i),temp_old(i),i)
  enddo
  return
end subroutine diffusion_check 
!--------1---------2---------3---------4---------5---------6---------7---------8
subroutine calc_diffusion_length(time,dt,temp,temp_old,ir) 
!--------1---------2---------3---------4---------5---------6---------7---------8
  use global_constants, only : temperature_Kelvin
  use global_constants, only : year, Myr
  use heat_radiogenic, only : iAl, ftime
  integer,intent(in) :: ir
  real(rkind),intent(in) :: time, dt, temp, temp_old
  real(rkind) :: iRT, diffusivity
  iRT = 1._rkind/(gas_constant*temp)
  diffusivity = calc_diffusivity(iRT)
  diffusion_length(ir) = diffusion_length(ir) + diffusivity*dt
  return
!--------1---------2---------3---------4---------5---------6---------7---------8
contains
real(rkind) function calc_diffusivity(iRT)
  real(rkind),parameter :: dfe1 = 10._rkind**(-8.34_rkind)
  real(rkind),parameter :: dfe2 = -338.e3_rkind
  real(rkind),intent(in) :: iRT 
  real(rkind) cal_diffusivity
  calc_diffusivity = dfe1*exp(dfe2*iRT)
  return 
end function calc_diffusivity
!--------1---------2---------3---------4---------5---------6---------7---------8
end subroutine calc_diffusion_length 
!--------1---------2---------3---------4---------5---------6---------7---------8
end module diffusion