!--------1---------2---------3---------4---------5---------6---------7---------8 module heat_parameters !--------1---------2---------3---------4---------5---------6---------7---------8 use global_constants, only : rkind use global_parameters, only : irmax implicit none private public :: set_surface_constant, coef_surface_constant, surface_constant & & , set_heat_parameters, thermal_conductivity, thermometric_diffusivity & & , specific_heat, latent_heat_init, latent_heat_vapor real(rkind) :: coef_surface_constant, surface_constant & & , thermal_conductivity(0:20), thermometric_diffusivity(0:20) & & , expansion_coefficient(0:20) , viscosity(0:20) , latent_heat(0:20) & & , specific_heat(0:20), latent_heat_init, latent_heat_vapor contains !--------1---------2---------3---------4---------5---------6---------7---------8 subroutine set_heat_parameters use global_constants, only : volume_ratio_rock, density, density_rock real(rkind),parameter :: thermal_conductivity_vapor = 1.e-2_rkind real(rkind),parameter :: thermal_conductivity_water = 0.56_rkind real(rkind),parameter :: thermal_conductivity_ice0 = 2.2_rkind real(rkind),parameter :: thermal_conductivity_ice = 3.69_rkind real(rkind),parameter :: thermal_conductivity_rock = 0.3248_rkind real(rkind),parameter :: specific_heat_vapor = 2000._rkind real(rkind),parameter :: specific_heat_water = 4200._rkind real(rkind),parameter :: specific_heat_ice0 = 1925._rkind real(rkind),parameter :: specific_heat_ice = 1380._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 :: specific_heat_rock = specific_heat_rockA real(rkind),parameter :: thermometric_diffusivity_rockA = & & thermal_conductivity_rock/(specific_heat_rock*density_rock) real(rkind),parameter :: thermometric_diffusivity_rockB = 0._rkind real(rkind),parameter :: expansion_coefficient_water = 2.1e-3_rkind real(rkind),parameter :: viscosity_water = 1.792e-6_rkind real(rkind),parameter :: latent_heat_ice = 3.34e5_rkind real(rkind),parameter :: latent_heat_water_vapor = 2.2e6_rkind real(rkind) :: mass_ratio_rock_mn6, mass_ratio_rock_mn15, mass_ratio_rock_mn16 thermal_conductivity(0) = thermal_conductivity_vapor thermal_conductivity(1) = thermal_conductivity_water thermal_conductivity(2) = thermal_conductivity_ice thermal_conductivity(4) = thermal_conductivity_rock thermal_conductivity(6) = thermal_conductivity_rock*volume_ratio_rock(6) & & + thermal_conductivity_ice*(1._rkind-volume_ratio_rock(6)) thermal_conductivity(15) = thermal_conductivity_rock*volume_ratio_rock(15) & & + thermal_conductivity_water*(1._rkind-volume_ratio_rock(15)) thermal_conductivity(16) = thermal_conductivity_rock*volume_ratio_rock(16) & & + thermal_conductivity_water*(1._rkind-volume_ratio_rock(16)) specific_heat(0) = specific_heat_vapor specific_heat(1) = specific_heat_water specific_heat(2) = specific_heat_ice specific_heat(4) = specific_heat_rock mass_ratio_rock_mn6 = volume_ratio_rock(6)*density(4)/density(6) specific_heat(6) = specific_heat_rock*mass_ratio_rock_mn6 & & + specific_heat_ice*(1._rkind-mass_ratio_rock_mn6) mass_ratio_rock_mn15 = volume_ratio_rock(15)*density(4)/density(15) specific_heat(15) = specific_heat_rock*mass_ratio_rock_mn15 & & + specific_heat_water*(1._rkind-mass_ratio_rock_mn15) mass_ratio_rock_mn16 = volume_ratio_rock(16)*density(4)/density(16) specific_heat(16) = specific_heat_rock*mass_ratio_rock_mn16 & & + specific_heat_water*(1._rkind-mass_ratio_rock_mn16) thermometric_diffusivity(1) = thermal_conductivity(1) & & /(density(1)*specific_heat(1)) thermometric_diffusivity(2) = thermal_conductivity(2) & & /(density(2)*specific_heat(2)) thermometric_diffusivity(4) = thermal_conductivity(4) & & /(density(4)*specific_heat(4)) thermometric_diffusivity(6) = thermal_conductivity(6) & & /(density(6)*specific_heat(6)) thermometric_diffusivity(15) = thermal_conductivity(15) & & /(density(15)*specific_heat(15)) thermometric_diffusivity(16) = thermal_conductivity(16) & & /(density(16)*specific_heat(16)) expansion_coefficient(1) = expansion_coefficient_water viscosity(1) = viscosity_water latent_heat(2) = latent_heat_ice latent_heat_init = latent_heat(2) latent_heat_vapor = latent_heat_water_vapor end subroutine set_heat_parameters !--------1---------2---------3---------4---------5---------6---------7---------8 subroutine set_surface_constant use global_constants, only : stefan_boltzman_constant, emissivity use global_parameters, only : temperature_initial, dr, iau surface_constant = 4._rkind*stefan_boltzman_constant*emissivity & & /thermal_conductivity(6)*temperature_initial(iau)**3._rkind coef_surface_constant = 2._rkind*surface_constant*dr end subroutine set_surface_constant !--------1---------2---------3---------4---------5---------6---------7---------8 end module heat_parameters !--------1---------2---------3---------4---------5---------6---------7---------8 module heat_radiogenic !--------1---------2---------3---------4---------5---------6---------7---------8 use global_constants, only : rkind, MeV_J, year implicit none private public :: heat_radiogenic_Al, iAlmax, heatAl, iAl, lambdaAl, initial_ratioAl & & , initial_ratioFe, initial_ratioK & & , ftime, mass_ratio_Al, mass_ratio_Al_hydrous & & , heat_radiogenic_Fe, heatFe, lambdaFe, mass_ratio_Fe, mass_ratio_Fe_hydrous & & , heat_radiogenic_K, heatK, lambdaK, mass_ratio_K, mass_ratio_K_hydrous integer,parameter :: iAlmax = 400 real(rkind),parameter :: halflife_Al = 7.2e5_rkind*year real(rkind),parameter :: standard_Al = 5.25e-5_rkind real(rkind),parameter :: mass_ratio_Al = 0.85e-2_rkind real(rkind),parameter :: mass_ratio_Al_hydrous = mass_ratio_Al/1.1_rkind real(rkind),parameter :: atomic_weight_Al = 26._rkind real(rkind),parameter :: decay_energy_Al = 3.16_rkind*MeV_J real(rkind),parameter :: halflife_Fe = 2.6e6_rkind*year real(rkind),parameter :: standard_Fe = 1.e-7_rkind real(rkind),parameter :: mass_ratio_Fe = 0._rkind real(rkind),parameter :: mass_ratio_Fe_hydrous = mass_ratio_Fe/1.1_rkind real(rkind),parameter :: atomic_weight_Fe = 60._rkind real(rkind),parameter :: decay_energy_Fe = 2.71_rkind*MeV_J real(rkind),parameter :: halflife_K = 1.248e9_rkind*year real(rkind),parameter :: standard_K = 1.559e-3_rkind real(rkind),parameter :: mass_ratio_K = 0._rkind real(rkind),parameter :: mass_ratio_K_hydrous = mass_ratio_K/1.1_rkind real(rkind),parameter :: atomic_weight_K = 40._rkind real(rkind),parameter :: decay_energy_K = 0.6550_rkind*MeV_J integer :: iAl real(rkind) :: heatAl(1:iAlmax), lambdaAl, initial_ratioAl(1:iAlmax) & & , initial_ratioFe(1:iAlmax), initial_ratioK(1:iAlmax) & & , ftime(1:iAlmax), heatFe(1:iAlmax), lambdaFe, heatK(1:iAlmax), lambdaK contains subroutine heat_radiogenic_Al use global_constants, only : year, mass_Hatom implicit none real(rkind),parameter :: latest = 5.e-7_rkind real(rkind) :: dlogAl, logAl, coef_decay dlogAl = (log10(standard_Al) - log10(latest))/real(iAlmax) lambdaAl = log(2._rkind)/halflife_Al coef_decay = decay_energy_Al*lambdaAl/(atomic_weight_Al*mass_Hatom) do iAl = 1, iAlmax logAl = log10(latest) + dlogAl*real(iAl) initial_ratioAl(iAl) = 10._rkind**logAl heatAl(iAl) = coef_decay*initial_ratioAl(iAl) ftime(iAl) = halflife_Al/log(2._rkind)*log(standard_Al/initial_ratioAl(iAl))/year enddo return end subroutine heat_radiogenic_Al subroutine heat_radiogenic_Fe use global_constants, only : year, mass_Hatom implicit none real(rkind) :: dlogAl, logAl, coef_decay lambdaFe = log(2._rkind)/halflife_Fe coef_decay = decay_energy_Fe*lambdaFe/(atomic_weight_Fe*mass_Hatom) do iAl = 1, iAlmax initial_ratioFe(iAl) = standard_Fe*exp(-lambdaFe*ftime(iAl)*year) heatFe(iAl) = coef_decay*initial_ratioFe(iAl) enddo return end subroutine heat_radiogenic_Fe subroutine heat_radiogenic_K use global_constants, only : year, mass_Hatom implicit none real(rkind) :: dlogAl, logAl, coef_decay lambdaK = log(2._rkind)/halflife_K coef_decay = decay_energy_K*lambdaK/(atomic_weight_K*mass_Hatom) do iAl = 1, iAlmax initial_ratioK(iAl) = standard_K*exp(-lambdaK*ftime(iAl)*year) heatK(iAl) = coef_decay*initial_ratioK(iAl) enddo return end subroutine heat_radiogenic_K end module heat_radiogenic